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Abstract 

Speeding up the data acquisition is one of the central aims to advance tomographic 
imaging. On the one hand, this reduces motion artifacts due to undesired movements, 
and on the other hand this decreases the examination time for the patient. In this article, 
we propose a new scheme for speeding up the data collection process in photoacoustic 
tomography. Our proposal is based on compressed sensing and reduces acquisition time 
and system costs while maintaining image quality. As measurement data we use random 
combinations of pressure values that we use to recover a complete set of pressure data 
prior to the actual image reconstruction. We obtain theoretical recovery guarantees 
for our compressed sensing scheme and support the theory by reconstruction results on 
simulated data as well as on experimental data. 

Keywords. Photoacoustic imaging, computed tomography, compressed sensing, lossless 
expanders, wave equation. 
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1 Introduction 


Photoacoustic tomography (PAT) is a recently developed non-invasive medical imaging tech¬ 
nology whose benefits combine the high contrast of pure optical imaging with the high spatial 
resolution of pure ultrasound imaging Hisiiisa. In order to speed up the measurement 
process, in this paper we propose a novel compressed sensing approach for PAT that uses 
random combinations of the induced pressure as measurement data. The proposed strategy 
yields recovery guarantees and furthermore comes with an efficient numerical implementation 
allowing high resolution real time imaging. We thereby focus on a variant of PAT using in¬ 
tegrating line detectors proposed in [T0l[53l[6]. Our strategy, however, can easily be adapted 
to more classical PAT setups using arrays of point-like detectors. 

Our proposal is based on the main components of compressed sensing, namely randomness 
and sparsity. Compressed sensing is one of the most influential discoveries in applied math¬ 
ematics and signal processing of the past decade [151 HI] • By combining the benefits of data 
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Figure 1.1: Basic principle of PAT. A semi-transparent sample is illuminated with a 
short optical pulse that induces an acoustic pressure wave. The induced pressure is measured 
outside of the sample and used to recover an image of the interior. 


compression and data acquisition it allows to recover a signal from far fewer linear measure¬ 
ments than suggested by Shannon’s sampling theorem. It has led to several new proposed 
sampling strategies in medical imaging, for example for speeding up MRI data acquisition 
(see SZlllH]) or completing under-sampled CT images m- Another prominent application 
of compressed sensing is the single pixel camera (see (22] ) that circumvents the use of several 
expensive high resolution sensors in digital photography. 


1.1 Photoacoustic tomography (PAT) 

PAT is based on the generation of acoustic waves by illuminating a semi-transparent sample 
with short optical pulses (see Figure ITT]) . When the sample is illuminated with a short laser 
pulse, parts of the optical energy become absorbed. Due to thermal expansion a subsequent 
pressure wave is generated depending on the structure of the sample. The induced pressure 
waves are recorded outside of the sample and used to recover an image of the interior, see 

[3 SSI ESI EH- 



Figure 1.2: PAT with integrating 
LINE DETEGTORS. An array of integrat¬ 
ing line detectors measures integrals of the 
induced acoustic pressure over a certain 
number of parallel lines. These data are 
used to recover a linear projection of the 
object in the first step. By rotating the ar¬ 
ray of line detectors around a single axis 
several projection images are obtained and 
used to recover the actual three dimen¬ 
sional object in a second step. 


In this paper, we consider a special variant of photoacoustic tomography that uses integrating 
line detectors for recording the pressure waves, as proposed in m- As illustrated in Figure 
o an array of line detectors is arranged around the investigated sample and measures 
integrals of the pressure wave over a certain number of parallel lines. Assuming constant 
speed of sound, the pressure integrated along the direction of the line detectors satisfies the 
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two dimensional wave equation 


{ d^p{x, t) — Ap{x, t) = 0 , for (x, <) S x (0, oo) 

p{x,0) = f{x), for X G (1) 

dtp (x, 0) = 0 , for X G , 

where the time scaling is chosen in such a way that the speed of sound is normalized to 
one. The initial datum / in m is the two dimensional projection image of the actual, three 
dimensional initial pressure distribution. 

Image reconstruction in PAT with integrating line detectors can be performed via a two-stage 
approach PEg. In the first step the measured pressure values corresponding to values of the 
solution of (HD outside the support of / are used to reconstruct a linear projection (the initial 
data in ([1])) of the three dimensional initial pressure distribution. This procedure is repeated 
by rotating the array of line detectors around a single axis which yields projection images 
for several angles. In a second step these projection images are used to recover the three 
dimensional initial pressure distribution by inverting the classical Radon transform. In this 
work we focus on the first problem of reconstructing the (two dimensional) initial pressure 
distribution / in (HD- 

Suppose that the integrating line detectors are arranged on the surface of a circular cylinder 
of radius p > 0, and that the object is located inside that cylinder (see Figure [T^)- The data 
measured by the array of line detectors are then modeled by 


Pj := 


p{zj, •) : [0,2p] -)> R, 

/pcos (27r(j - l)/N)\ 
VPsin(27r(j - 1)/^)/ ’ 


for j = 1,..., , 


( 2 ) 

(3) 


where pj is the pressure signal corresponding to the j-th line detector. Since the two dimen¬ 
sional initial pressure distribution / is supported in a disc of radius p and the speed of sound 
is constant and normalized to one, no additional information is contained in data p {zj , t) for 
t > 2p. This can be seen, for example by exploiting the explicit relations between the two 
dimensional pressure signals p{zj, ■) and the spherical means (compare with Subsection l3.3l) 
of the initial pressure distribution; see M- 

Of course, in practice also the functions pj: [0, 2p] —>• R have to be represented by discrete 
samples. However, temporal samples can easily be collected at a high sampling rate compared 
to the spatial sampling, where each sample requires a separate sensor. It is therefore natural 
to consider the semi-discrete data model (ED- Our compressed PAT scheme could easily be 
adapted to a fully discretized data model. 


1.2 Compressed sensing PAT 

When using data of the form ED, high resolution imaging requires the number N of detector 
locations to be sufficiently large. As the fabrication of an array of parallel line detectors is 
demanding, most experiments using integrating line detectors have been carried out using a 
single line detector, scanned on circular paths using scanning stages EH ED]. Recently, systems 
using arrays of 64 parallel line detectors have been demonstrated [MIE]. The most costly 
building blocks in such devices are the analog to digital converters (ADC). For completely 
parallel readout, a separate ADC is required for every detector channel. In order to reduce 
costs, in these practical implementations two to four line detector channels are multiplexed to 
one ADC. For collecting the complete pressure data, the measurements have to be performed 
two (respectively four) times, because only 32 (respectively 16) of the 64 line measurements 
can be read out in parallel. This, again, leads to an increased overall measurement time. For 
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example, using an excitation laser with a repetition rate of 10 Hz and two times multiplexing, 
a measurement time of 0.2 s for a projection image has been reported in m- Without 
multiplexing, this measurement time would reduce to the half. 

In order to speed up the scanning process and to reduce system costs, in this paper we propose 
a novel compressed sensing approach that allows to perform a smaller number of random 
measurements with a reduced number of ADCs, while retaining high spatial resolution. For 
that purpose, instead of collecting individually sampled data pj (t) as in ([2]) , we use random 
combinations 

».(t) = E Pj (t) for i £ {1,..., m} and t £ [0, 2p] , (4) 

j&Ji 

where m <C is the number of compressed sensing measurements and Ji C {1,...,A^} 
corresponds to the random set of detector locations contributing to the i-th measurement. 
In the reconstruction process the random linear combinations yi{t) are used to recover the 
full set of pressure data pi{t),... ,pN{t) using compressed sensing techniques. The initial 
pressure distribution / in © is subsequently recovered from the completed pressure data by 
applying standard PAT reconstruction algorithms such as time reversal [miMiiiD] or filtered 
backprojection [SJ |221 HHIHHIIHIS2] ■ 

A naive approach for recovering the pressure data from the random measurements yi(t) 
would be to solve for Pj{t) separately for each t £ [0,2p]. Since m N, this is a 
severely underdetermined system of linear equations and its solution requires appropriate 
prior knowledge of the unknown parameter vector. Compressed sensing suggests to use the 
sparsity of the parameter vector in a suitable basis for that purpose. However, recovery 
guarantees for zero/one measurements of the type (|3]) are basis-dependent and require the 
parameter to be sparse in the standard basis rather than sparsity in a different basis such 
as orthonormal wavelets (see Subsection 12.21) . However, for pressure signals (H)) of practical 
relevance such sparsity assumption in the original basis does not hold. 

In this work we therefore propose a different approach for solving by exploiting special 
properties of the data in PAT. For that purpose we apply a transformation that acts in the 
temporal variable only, and makes the transformed pressure values sufficiently sparse in the 
spatial (angular) component. In Subsection Id.dl we present an example of such a transform. 
The application of a sparsifying transform to ® yields linear equations with unknowns being 
sparse in the angular variable. It therefore allows to apply sparse recovery results for the 
zero/one measurements under consideration. 

1.3 Relations to previous work 

A different compressed sensing approach for PAT has been considered in [55l [31] . In these 
articles, standard point samples (such as ([H) have been used as measurement data and no 
recovery guarantees have been derived. Further, in [551 [5T] the phantom is directly recon¬ 
structed from the incomplete data, whereas we first complete the data using sparse recovery 
techniques. Our approach is more related to a compressed sensing approach for PAT using 
a planar detector array that has been proposed in [35] and also uses random zero/one com¬ 
binations of pressure values and recovers the complete pressure prior to the actual image 
reconstruction. However, in |35] the sparsifying transform is applied in spatial domain where 
recovery guarantees are not available as noted above. We finally notice that our proposal of 
using a sparsifying temporal transform can easily be extended to planar detector arrays in 
two or three spatial dimensions; compare Section [5] 
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1.4 Outline 


The rest of this paper is organized as follows. In Section [5] we review basic results from 
compressed sensing that we require for our proposal. We therefore focus on recovery guar¬ 
antees for zero/one matrices modeled by lossless expanders; see Subsection 12.21 In Section |3] 
we present the mathematical framework of the proposed PAT compressed sensing scheme. 
The sparsity in the spatial variable, required for ^^-minimization, is obtained by applying 
a transformation acting in the temporal variable. An example of such a transformation is 
given in Subsection 13.31 In Section |T] we present numerical results supporting our theoretical 
investigations. The paper concludes with a short discussion in Section [SI 


2 Background from compressed sensing 

In this section we shortly review basic concepts and results of compressed sensing (sometimes 
also termed compressive sampling). Our main focus will be on recovery results for lossless 
expanders, which are the basis of our PAT compressed sensing scheme. 

2.1 Compressed sensing 

Suppose one wants to sense a high dimensional data vector x = {xi,... ,xn) G K.'^, such 
as a digital image. The classical sampling approach is to measure each component of Xi 
individually. Hence, in order to collect the whole data vector one has to perform N separate 
measurements, which may be too costly. On the other hand it is well known and the basis of 
data compression algorithms, that many datasets are compressible in a suitable basis. That 
is, a limited amount of information is sufficient to capture the high dimensional vector x. 

Compressed sensing incorporates this compressibility observation into the sensing mechanism 
[13 nil m]. Instead of measuring each coefficient of the data vector individually, one collects 
linear measurements 

Ax = y, (5) 

where A S is the measurement matrix with m N, and y = {yi,... ,ym) € K.™ is 

the measurement vector. Any component of the data vector can be interpreted as a scalar 
linear measurement performed on the unknown x, and the assumption m <C fV means that 
far fewer measurements than parameters are available. As to <C fV, the system ((3) is highly 
underdetermined and cannot be uniquely solved (at least without additional information) by 
standard linear algebra. 

Compressed sensing overcomes this obstacle by utilizing randomness and sparsity. Recall that 
the vector x = (xi ,... ,xn) is called s-sparse if the support 

supp (x) := {j G {1,...,N} : Xj ^ 0} (6) 

contains at most s elements. Results from compressed sensing state that for suitable A, any 
s-sparse x G R" can be found via the optimization problem 

N 

minimize llzIL = > \zA 

zeR'^ " (7) 

such that Az = y . 

By relaxing the equality constraint Az = y, the optimization problem ([7]) can be adapted to 
data which are only approximately sparse and noisy m- 
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A sufficient condition to guarantee recovery is the so called restricted isometry property (RIP), 
requiring that for any s-sparse vector x, we have 

(1 — (5)||x||2 < ||Ax ||2 < (1 + (5)||x|| 2 for some small 5 G (0,1). 

The smallest constant <5 satisfying this inequality is called the s-restricted isometry constant 
of A and denoted by 5s- Under certain conditions on 5s, recovery guarantees for sparse and 
approximately sparse data can be obtained, see for example [uni]- 

While the restricted isometry itself is deterministic, to date all constructions that yield near- 
optimal embedding dimensions m are based on random matrices. Sub-gaussian random 
matrices satisfy the RIP with high probability for an order-optimal embedding dimension 
m = 0{s\og{N/s)), see e.g. [T]. Partial random Fourier matrices (motivated by MRI mea¬ 
surements) and subsampled random convolutions (motivated by remote sensing) have been 
shown to allow for order-optimal embedding dimensions up to logarithmic factors, see |581I56] 
and [^SO], respectively. 

The sparsity is often not present in the standard basis of , but in a special sparsifying 
basis, such as wavelets. For matrices with subgaussian rows this does not cause a problem, 
as the rotation invariance of subgaussian vectors ensures that after incorporating an orthog¬ 
onal transform, the resulting random matrix construction still yields RIP matrices with high 
probability. As a consequence, the sparsifying basis need not be known for designing the 
measurement matrix A, which is often referred to as universality of such measurements [T]. 

Many structured random measurement systems including the partial random Fourier and 
the subsampled random convolution scenarios mentioned above, however, do not exhibit 
universality. For example, one can easily see that subsampled Fourier measurements cannot 
suffice if the signal is sparse in the Fourier basis. While it has been shown that this problem 
can be overcome by randomizing the column signs [41) , such an alteration often cannot 
be implemented in the sensing setup. Another way to address this issue is by requiring 
incoherence between the measurement basis and the sparsity basis m- That is, one needs 
that inner products between vectors of the two bases are uniformly small. If not all, but 
most of these inner products are small, one can still recover, provided that one adjusts the 
sampling distribution accordingly; this scenario includes the case of Fourier measurements 
and Haar wavelet sparse signals |42) . 

Incoherence is also the key to recovery guarantees for gradient sparse signals. Namely, many 
natural images are observed to have an approximately sparse discrete gradient. As a conse¬ 
quence, it has been argued using a commutation argument that one can recover the signal 
from uniformly subsampled Fourier measurements via minimizing the norm of the discrete 
gradient, the so-called total variation (TV) [TS]. TV minimization had already proven to 
be a favorable method in image processing, see, for example HH]. A problem with this 
approach is that the compressed sensing recovery guarantees then imply good recovery of 
the gradient, not of the signal itself. Small errors in the gradient, however, can correspond 
to substantial errors in the signal, which is why this approach can only work if no noise is 
present. A refined analysis that allows for noisy measurements requires the incoherence of 
the measurement basis to the Haar wavelet basis m- Again, TV minimization is consid¬ 
ered for recovery. By adjusting the sampling distribution, these results have also been shown 
to extend to Fourier measurements and other systems with only most measurement vectors 
incoherent to the Haar basis [42]. 

For the measurement matrices considered in this work, namely zero/one matrices based on ex¬ 
pander graphs, recovery guarantees build on an .^^-version of the restricted isometry property, 
namely one requires that 

(I-5)||x||i<||Ax||i<(l + 5)||x||i 

for all sufficiently sparse x and some constant i5 > 0; see [7] and Subsection 12.21 below. As the 
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£^-norm is not rotation invariant, basis transformations typically destroy this property. That 
is, not even incoherence based recovery guarantees are available; recovery results only hold 
in the standard basis. Thus an important aspect of our work will be to ensure (approximate) 
sparsity in the standard basis. This will be achieved by applying a particular transformation 
in the time variable. 


2.2 Recovery results for lossless expanders 

Recall that we seek recovery guarantees for a measurement setup, where each detector is 
switched on exactly d out of m times. That is, one obtains a binary measurement matrix A G 
{0,1}™^^ with exactly d ones in each column. It therefore can be interpreted as the adjacency 
matrix of a left d-regular bipartite graph. Under certain additional conditions, such a bipartite 
graph is a lossless expander (see Definition 12.1|) which, as we will see, guarantees stable 
recovery of sparse vectors. Expander graphs have been used since the 1970s in theoretical 
computer science, originating in switching theory from modeling networks connecting many 
users, cf. |39j for further applications. They have also been useful in measure theory, where 
it was possible to solve the Banach-Ruziewicz problem using tools from the construction of 
expander graphs, see |45] for a detailed examination of this connection. For a survey on 
expander graphs and their applications, see for example |M1146) . 

Compressed sensing with expander graphs has been considered in [31 [Ml [SHIM], where also 
several efficient algorithms for the solution of compressed sensing problems using expander 
graphs have been proposed. A short review of sparse recovery algorithms using expander 
like matrices is given in |26j . In this subsection we recall main compressed sensing results 
using lossless expanders as presented in the recent monograph [25) . where the proofs of all 
mentioned theorems can be found in Section 13. 


N 



Figure 2.1: Example of a left d-regular bipartite graph with d = 2, N = 7 left vertices, and 
TO = 5 right vertices. Here d = 2 because exactly 2 edges emerge at each left vertex. For 
J = {2, 3} the set of right vertices connected to J is given by R{J) = {2, 3,4}. 

Recall that a bipartite graph consists of a triple (L, R, E), where L is the set of left vertices, R 
the set of right vertices, and E C L x R is the set of edges. Any element (j, i) £ E represents 
an edge with left vertex j € L and right vertex i £ R. 

A bipartite graph {L,R,E) is called d-regular for some d > 0, if for every given left vertex 
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j G L, the number of edges (j, i) € E emerging from j is exactly equal to d. Finally, for any 
subset J C L of left vertices, let 

R{J) := {i G R: there exists some j G J with (j,i) G E} 

denote the set of right vertices connected to J. Obviously, for any left d-regular bipartite 
graph and any J C L, we have |i?(J)| < d| J|- 

Definition 2.1 (Lossless expander). Let {L,R,E) be a left d-regular bipartite graph, s G N, 
and 6 G [0,1]. Then (L, R, E) is called an (s, d, 0)-lossless expander, if 

|^(•^)| > (1 ~ 0)d\J\ for all J C L with | J| < s . (8) 

For any left d-regular bipartite graph, the smallest number 6 > 0 satisfying is called the 
s-th restricted expansion constant and denoted by 9s. 


The following theorem states that a randomly chosen d-regular bipartite graph will be a 
lossless expander with high probability. 


Theorem 2.2 (Regular bipartite graphs are expanders with high probability). For every 
0 < e < every 9 G (0,1) and every s G N, the proportion of {s,d, 9)-lossless expanders 
among the set of all left d-regular bipartite graphs having N left vertices and m right vertices 
exceeds 1 — e, provided that 


d = 



and 


m > cgs In 



(9) 


Here eg is a constant only depending on 9, e is Euler’s constant, ln( •) denotes the natural 
logarithm and [a;] denotes the smallest integer larger or equal to x. 


According to Theorem 12.21 any randomly chosen left d-regular bipartite graph is a lossless 
expander with high probability, provided that (|9]) is satisfied. The following theorem states 
that the adjacency matrix A G {0, ^ 

{j,i)GE, (10) 

of any lossless expander with left vertices L = {1,..., A} and right vertices i? = {1, ..., m} 
yields stable recovery of any sufficiently sparse parameter vector. The result was first estab¬ 
lished in [7], we present the version found in [25] . 

Theorem 2.3 (Recovery guarantee for lossless expanders). Let A G {0 ,be the adja¬ 
cency matrix of a left d-regular bipartite graph having 92s <1/6. Further, let x G C^, ij > 0 
and e G C™ satisfy ||e||i < rj, set b := Ax -|- e, and denote by x* a solution of 


minimize llzIL 

zGC« ^ 

such that II Az — b||^ < ? 7 . 


( 11 ) 


Then 


l|x 


X*||l < 


2 ( 1 - 202 .) 
(1 - 602 .) 


^.(x)! 


4 

(l- 602.)d'^' 


Here the quantity tTs(x)i := inf{||x — z||i: z is s-sparse} measures by how much the vector 
X G fails to he s-sparse. 


Combining Theorems l2.2l and l2.31 we can conclude that the adjacency matrix A of a randomly 
chosen left d-regular bipartite graph will, with high probability, recover any sufficiently sparse 
vector X G by basis pursuit reconstruction (HH). 
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3 Mathematical framework of the proposed PAT com¬ 
pressed sensing scheme 

In this section we describe our proposed compressed sensing strategy. As mentioned in the 
introduction, we focus on PAT with integrating line detectors, which is governed by the two 
dimensional wave equation (ED- In the following we first describe the compressed sensing mea¬ 
surement setup in Subsection 111. II and describe the sparse recovery strategy in Subsection ld.2l 
As the used pressure data are not sparse in the original domain we introduce a temporal 
transform that makes the data sparse in the spatial domain. In Subsection 13.31 we present an 
example of such a sparsifying temporal transform. In Subsection 13.41 we finally summarize 
the whole PAT compressed sensing scheme. 


Figure 3.1: Left: Classical PAT sampling, where a single detector is moved around the 
object to collect individual pressure signal pj. Right: Compressed sensing approach where 
each measurement consists of a random zero/one combination of individual pressure values. 




3.1 Compressed sensing PAT 

We define the unknown full sample pressure pj: [0, 2p] —R, for j S {1, ..., N} by (I2|), ([3]), 
where p is the solution of the two dimensional wave equation (ED- We suppose that the (two 
dimensional) initial pressure distribution / is smooth and compactly supported in i3_R(0) 
which implies that any pj is smooth and vanishes in a neighborhood of zero. Furthermore we 
assume that the data ([2]) are sampled finely enough to allow for / to be reconstructed from 
Pj by standard PAT reconstruction algorithms such as time reversal [TTl [331IBO] or filtered 
backprojection [SJ [331 [331 [3311331IM] ■ 

Instead of measuring each pressure signal pj separately, we take m compressed sensing mea¬ 
surements. For each measurement, we select sensor locations Zj with j S Ji at random and 
take the sum of the corresponding pressure values pj. Thus the i-th measurement is given by 

■=^pA'^) for f G {1,... ,m} , (12) 

j&Ji 

where Ji C {1,... ,fV} corresponds to the set of all detector locations selected for the f-th 
measurement, and m ^ N is the number of compressed sensing measurements. 

In practice, the compressed sensing measurements could be realized by summation of several 
detector channels, using a configurable matrix switch and a summing amplifier. Even more 
simply, a summing amplifier summing over all N channels could be used, while individual 
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detector channels Zj are turned on or off, using solid state relays. Thereby, only one ADC 
is required for one compressed sensing measurement. Performing m compressed sensing mea¬ 
surements in parallel is facilitated by using m ADCs in parallel, and the according number 
of matrix switches, relays, and summing amplifiers. 

In the following we write 

p: (0,oo)-)> t p(t) = (pi(t),... ,pAr(t))'^ , (13) 

y: (0,oo)R™: t y(t) = ( 2 /i(t),... , (14) 


for the vector of unknown complete pressure signals and the vector of the corresponding 
measurement data, respectively. Further, we denote by 

A £ {0, 1 }™^^ with entries := i J ^ * 

' \0, forj^J,, 

the matrix whose entries in the i-th row correspond the sensor locations selected for the i-th 
measurement. In order to apply the exact recovery guarantees from Subsection l2.2l we require 
that each row of A contains exactly d ones, where d £ N is some fixed number. Practically, 
this means that each detector location contributes to exactly d of the m measurements, which 
also guarantees that the measurements are well calibrated and there is no bias towards some 
of the detector locations. 

Recovering the complete pressure data @ from the compressed sensing measurements ® 
can be written as an uncoupled system of under-determined linear equations. 


Ap(t) = y{t) for any t £ [0, 2p ], (15) 

where Ap(t) := A (p(t)) for any t. From dTSl) . we would like to recover the complete set of 
pressure values p{t) for all times t £ [0, 2p]. Compressed sensing results predict that under 
certain assumptions on the matrix A any s-sparse vector p(t) can be recovered from Ap(t) 
by means of sparse recovery algorithms like ^^-minimization. 

Similar to many other applications (cf. Section [5] above), however, we cannot expect sparsity 
in the original domain. Instead, one has p(t) = 'l'x(t), where 4/ is an appropriate orthogonal 
transform and x(t) £ R'^ is a sparse coefficient vector. This yields a sparse recovery problem 
for x(t) involving the matrix Ad^, which does not inherit the the recovery guarantees of 
Subsection 12.21 Hence we have to find a means to establish sparsity without considering 
different bases. Our approach will consist of applying a transformation in the time domain 
that sparsifies the pressure in the spatial domain. A further advantage of working in the 
original domain is that the structure of A allows the use of specific efficient algorithms like 
sparse matching pursuit [5] or certain sublinear-time algorithms like [371 Ell- 


3.2 Reconstruction strategy 

We denote by t/([0, 2p]) the set of all infinitely differentiable functions g: [0,2p] —?> R that 
vanish in a neighborhood of zero. To obtain the required sparsity, we will work with a 
sparsifying transformation 

T: e([0,2p])^5([0,2p]), (16) 

that is, Tp(t) £ R'^ can be sufficiently well approximated by a sparse vector for any t £ [0, 2p] 
and certain classes of practically relevant data p{t). Here we use the convention that T applied 
to a vector valued function g = {gi,..., gk) is understood to be applied in each component 
separately, that is. 


Tg(t):=((Tgi)(t),...,(Tgfc)(f)) , fort£[0,2p]. (17) 
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We further require that T is an injective mapping, such that any g G 0{[0,2p]) can be 
uniquely recovered from the transformed data Tg. See Subsection Id.31 for the design of such 
a sparsifying transformation. 

Since any temporal transformation interchanges with A, application of the sparsifying tem¬ 
poral transformation T to the original system (ESI yields 

A (Tp(t)) = Ty(t) for any t G [0, 2p]. (18) 

According to the choice of the temporal transform, the transformed pressure Tp(t) can be 
well approximated by a vector that is sparse in the spatial component. We therefore solve, 
for any t G [OiT], the following ^^-minimization problem 


minimize ||q||, 

qeK" (19) 

such that ||Aq-Ty(t)||,<, 7 , 


for some error threshold rj > 0. As follows from Theorem [231 the solution q*(t) of the 
minimization problem m provides an approximation to Tp(t) and consequently we have 
T~^q*(t) ~ p(t) for all t G [0,2/9]. Note that the use of ^^-norm || • ||^ for the constraint in 
(fT9l is required for the stable recovery guarantees for our particular choice of the matrix A 
containing zero/one entries using data that is noisy and only approximately sparse. 

As a further beneht, compressed sensing measurements may show an increased signal to noise 
ratio. For that purpose consider Gaussian noise in the pressure data, where instead of the 
exact pressure data p(t), we have noisy data p(t) = p(t) -I- r). For the sake of simplicity 
assume that the entries pj ^ A/’(0,ct^) of rj are independent and identically distributed. The 
corresponding noisy (and rescaled) measurements are then given by 


1 

W\ 


E 


{Pj{t) + Pj) 


1 

W\ 


jeJi 


1 

W\ 


E^j- 


The variance in the compressed sensing measurements is therefore cr^/j J/j compared to in 
the individual data pj (t). Assuming some coherent averaging in the signal part this yields 
an increased signal to noise ratio reflecting the inherent averaging of compressed sensing 
measurements. 


3.3 Example of a sparsifying temporal transform 

As we have seen above, in order to obtain recovery guarantees for our proposed compressed 
scheme, we require a temporal transformation that sparsifies the pressure signals in the an¬ 
gular component. In this section, we construct an example of such a sparsifying transform. 

Since the solution of the two dimensional wave equation can be reduced to the spherical 
means, we will construct such a sparsifying transform for the spherical means 

M /(z, r) := ^ [ f{z + ruj)da{uj) for (z, r) € dBniO) x (0, oo ). 

In fact, the solution of the two dimensional wave equation El can be expressed in terms of 
the spherical means via p{z,t) = dt fg rM.f(z,r)/\/t^ — r^dr, see [35]. By using standard 
tools for solving Abel type equations, the last expression can be explicitly inverted, resulting 
in (Mf)(z,r) = 2/7r p(z, t)l\/r‘^ — (see ElllZ])- Hence any sparsifying transformation 
for the spherical means M / also yields a sparsifying transformation for the solution of the 
wave equation El, and vice versa. 
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We found empirically, that drV H,. 5^ M / is sparse in the spatial direction for any function / 
that is the superposition of few indicator functions of regular domains. Here drg denotes the 
derivative in the radial variable, 

g{z, r) = - [ , for (z, r) G x (0, oo), 

Jr r-s 

the Hilbert transform of the function g: 95^(0) x R —>■ R in the second component, and r 
the multiplication operator that maps the function (z,r) i—^ g{z,r) to the function (z,r) i—> 
rg{z, r). Further, the spherical means M / : dBfi{0) x R —?> R are extended to an odd function 
in the second variable. For a simple radially symmetric phantom the sparsity of drr M / 
is illustrated in Figure 13.21 Thus we can choose T = drV dr as a sparsifying temporal 
transform for the spherical means. 



Figure 3.2: Sparsity induced by dr{rllr dr ). The left image shows the simple disc-like 
phantom / (characteristic function of a disc), the middle image shows the filtered spherical 
means r 9^ M / and the right image shows the sparse data dr {r Hr 9^ M /). 


The function r 9r M / also appears in the following formula for recovering a function from 
its spherical means derived in [23) . 

Theorem 3.1 (Exact reconstruction formula for spherical means). Suppose f G C°“(R^) is 
supported in the closure o/H/j(0). Then, for any x G we have 

/ (rHr9rM/) (z, |a; - z|)ds(z). (20) 

27rb: JdBuio) 

Proof. See [531 Corollary 1.2]. □ 

In the practical implementation the spherical means are given only for a discrete number of 
centers Zj G 9Hi^(0) yielding semi-discrete data similar to ([5]), ([3]). Formula (150]) can easily be 
adapted to discrete or semi-discrete data yielding a filtered backprojection type reconstruction 
algorithm; compare El Section 4]. So if we can find the filtered spherical means 

(r Hr dr M/) {zj, ■) for all detector locations Zj G 95^(0) (21) 

we can obtain the desired reconstruction of / by applying the backprojection operator (the 
outer integration in dini)) to rHr9rM/. 

3.4 Summary of the reconstruction procedure 

In this section, we combine and summarize the compressed sensing scheme for photoacoustic 
tomography as described in the previous sections. Our proposed compressed sensing and 
sparse recovery strategy takes the following form. 
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(csl) Create a matrix A G {0,as the adjacency matrix of a randomly selected left 
d-regular bipartite graph. That is, A is a random matrix consisting of zeros and ones 
only, with exactly d ones in each column. 

(cs2) Perform m measurements, whereby in the i-th measurement pressure signals corre¬ 
sponding to the nonzero entries in i-th row of A are summed up, see Equations o 
and m- This results in measurement data Ap(t) = y(t) G R™ for any t G [0, 2p]. 

(cs3) Choose a transform T acting in the temporal direction, which sparsifies the pressure 
data p along the spatial direction; compare Equation dm. 

(cs4) For any t G [0, 2p] and some given threshold rj, perform t'^-minimization (fT^ resulting 
in a sparse vector q*(t) satisfying ||Aq*(t) — Ty(t)||i < p. 

(cs5) Use p*(t) = T~^q*(t) as the input for a standard PAT inversion algorithm for complete 
data, such as time reversal or filtered backprojection. 


As we have seen in Subsection o the procedure |(csl)ff(cs5)| yields a close approximation 
to the original function / if the transformed data Tp(t) are sufficiently sparse in the spatial 
direction. The required sparsity level is hereby given by the expander-properties of the matrix 
A. Note that for exact data and exactly sparse data, we can use the error threshold p = 0. 
In the more realistic scenario of noisy data and Tp(t) being only approximately sparse, we 
solve the optimization problem m to yield a near optimal solution with error level bounded 
by the noise level. 


4 Numerical results 

To support the theoretical examinations in the previous sections, in this section we present 
some simulations using the proposed compressed sensing method. We first present reconstruc¬ 
tion results using simulated data and then show reconstruction results using experimental 
data. 

4.1 Results for simulated data 

As in Subsection E31 we work with the equivalent notion of the spherical means instead of 
directly working with the solution of the wave equation ©• In this case the compressed 
sensing measurements provide data 

yi{r) = rrij (r) for i G {I,..., m} and t G [0, 2p], (22) 

j&Ji 

where mj = {M.f){zj, •) denote the spherical means collected at the j-the detector loca¬ 
tion Zj. We further denote by m(t) = . .. ,mAr(t))^ the vectors of unknown complete 

spherical means and by y{t) = ..., ym{i)Y vector of compressed sensing measure¬ 
ment data. Finally, we denote by A G {0, the compressed sensing matrix such that 

(l2^ can be rewritten in the form Am = y. 

As proposed in Subsection l3.3l we use T = dr(r dr) as a sparsifying transform for the spher¬ 
ical means. An approximation to drV Hr 9r M / can be obtained from compressed sensing 
measurements in combination via .^^-minimization. For the recovery of the original function 
from the completed measurements, we use one of the inversion formulas of [23) presented 
given in Theorem 13.11 Recall that this inversion formulas can be implemented by applying 
the circular back-projection to the filtered spherical means r dr M /. 
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In order to obtain an approximation to the data m from the sparse intermediate reconstruc¬ 
tion SrCr Hf i9r M/)(•, r), one has to perform one numerical integration along the second 
dimension after the ^^-minimization process. We found that this numerical integration intro¬ 
duces artifacts in the reconstruction of ?’ 9,. M / , required for application of the inversion 

formula (EOl), see the middle image in Figure UTTJ These artifacts also yield some undesired 
blurring in the final reconstruction of /; see the middle image in Figure 



Figure 4.1: Reconstruction of filtered spherical means for N = 200 and m = 100. 
Left: Reconstruction of 9^1" Hf. Sr M / using £^-minimization. Center: Result of integrating 
the ^^-reconstruction in the radial direction. Right: Result of directly recovering r 9^ M / 
by TV-minimization. 

In order to overcome the artifact introduced by the numerical integration, instead of ap¬ 
plying an additional radial derivative to r dr to obtain sparsity in the spatial direc¬ 
tion, we will apply one-dimensional total variation minimization (TV) for directly recovering 
(r Hr 9r M /)(•, r). Thereby we avoid performing numerical integration on the reconstructed 
sparse data. Furthermore, this yields much better results in terms of image quality of the 
final image /. 

In our interpretation, this performance discrepancy is comparable to the difference between 
uniform and variable density samples for subsampled Fourier measurements. While |15j proves 
recovery of the discrete gradient, this does not carry over to the signal in a stable way - 
a refined analysis was required [sniiia. Similarly, we expect that a refined analysis to be 
provided in subsequent work can help explain the quality gap between and TV minimization 
that we observe in our scenario. 

In order to approximately recover r 9,. M / from the compressed sensing measurements, 
we perform, for any r G [0, 2p], one-dimensional discrete TV-minimization 

||Aq-(rHr9^y)(r)||2-f AIIqIItv^ mm . (23) 

Here ||q||xv = Iqj+i — Qj\ denotes the discrete total variation using the periodic 

extension qn+i ■= 9i- The one-dimensional total variation minimization problem (I23p can be 
efficiently solved using the fast iterative shrinkage thresholding algorithm (FISTA) of Beck and 
Teboulle The required proximal mapping for the total variation can be computed in 0{N) 
operation counts by the tautstring algorithm (see [iai2nii2H])- The approximate solution of 
(l23)l therefore only requires 0{NmNiter) floating point operations, with Nuer denoting the 
number of iterations in the FISTA. Assuming the radial variable to be discretized using 0{N) 
samples, the whole data completion procedure by (ESI only requires 0{N‘^mNiter) operation 
counts. Since we found that fewer than 100 iterations in the FISTA are often sufficient for 
accurate results, the numerical effort of data completion is only a few times higher than that 
of standard reconstruction algorithms in PAT. 

Figures UT] and EEl show results of simulation studies for a simple phantom, where the initial 
pressure distribution / is the characteristic function of a disc. For the compressed sensing 
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Figure 4.2: Reconstruction results for disc-like phantom using N = 200 and 
m = 100. Left: Reconstruction from 100 standard measurements. Center: Compressed sens¬ 
ing reconstruction using ^^-minimization. Right: Compressed sensing reconstruction using 
TV-minimization. The reconstruction from standard measurements contains under-sampling 
artifacts which are not present in the compressed sensing reconstructions. Further, the use 
of TV-minimization yields much less blurred results than the use of ^^-minimization. 


reconstruction, we used m = 100 random measurements instead of fV = 200 standard point 
measurements. As one can see from the right image in Figure 14.21 the completed data 
r Hr dr M /, for this simple phantom, is recovered almost perfectly from the compressed 
sensing measurements by means of TV-minimization. The reconstruction results in Figure|421 
show that the combination of our compressed sensing approach with TV-minimization yields 
much better results than the use of £^-minimization. Fur comparison purpose, the left image 
in Figure 14.21 shows the reconstruction from 100 standard measurements. Observe that the 
use of TO = 100 random measurements yields better result than the use of the 100 standard 
measurements, where artifacts due to spatial under-sampling are clearly visible. 


4.2 Results for real measurement data 

Experimental data have been acquired by scanning a single integrating line detector on a 
circular path around a phantom. A bristle with a diameter of 120 /im was formed to a 
knot and illuminated from two sides with pulses from a frequency doubled Nd:YAG laser 
with a wavelength of 532nm. The radiant exposure for each side was below 7.5™J/cm^. 
Generated photoacoustic signals have been detected by a graded index polymer optical fiber 
being part of a Mach-Zehnder interferometer, as described in [30] ■ Ultrasonic signals have 
been demodulated using a self-made balanced photodetector, the high-frequency output of 
which was sampled with a 12-bit data acquisition device. A detailed explanation of the used 
photo-detector and electronics can be found in [3]. The polymer fiber detector has been 
scanned around the phantom on a circular path with a radius of 6 mm and photoacoustic 
signals have been acquired on 121 positions. The scanning curve was not closed, but had an 
opening angle of 7r/2rad. Hence photoacoustic signals have been acquired between 7r/8rad 
and 157r/8rad. 

Using these experimental data, compressed sensing data have been generated, where each 
detector location was used d = 10 times and to = 60 measurements are made in total. The 
reconstruction of the complete measurement data has been obtained by one-dimensional dis¬ 
crete TV-minimization ([33| as suggested in the Section I3TT1 The measured and the recovered 
complete pressure data are shown in the top row in Figure 14.31 The bottom row in Fig¬ 
ure |T]3] shows the reconstruction results from 121 standard measurements (bottom left) and 
the reconstruction from 60 compressed sensing measurements (bottom center). Observe that 
there is only a small difference between the reconstruction results. This clearly demonstrates 
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the potential of our compressed sensing scheme (csl) - (cs5) for decreasing the number of 
measurements while keeping image quality. For comparison purpose we also display the re¬ 
construction using 60 standard measurement (bottom right). Compared to the compressed 
sensing reconstruction using the same number of measurements the use of standard measure¬ 
ments shows significantly more artifacts which are due to spatial under-sampling. 


X 10“® 



X 10"® 



0 ^ ' ' ' - 

-3-2-1 0 1 2 3 



X [mm] X [mm] x [mm] 

Figure 4.3: Reconstruction results for PAT measurements. Top left: Measured 
pressure data for N = 121 detector positions. Top right: Compressed sensing reconstruction 
of the full pressure data from m = 60 compressed measurements. Bottom left: Reconstruction 
from the full measurement data using N = 121 detector positions. Bottom center: Recon¬ 
struction from m = 60 compressed sensing measurements. Bottom right: Reconstruction 
from half of the measurement data using 60 detector positions. 


5 Discussion and outlook 

In this paper, we proposed a novel approach to compressive sampling for photoacoustic to¬ 
mography using integrating line detectors providing recovery guarantees for suitable datasets. 
Instead of measuring pressure data pj at any of the N individual line detectors, our approach 
uses m random combinations of pj with m <C as measurement data. The reconstruc- 


16 




















tion strategy consists of first recovering the pressure values pj for j G N} and then 

applying a standard PAT reconstruction for obtaining the final photoacoustic image. For re¬ 
covering the individual pressure data, we propose to apply a sparsifying transformation that 
acts in the temporal variable and makes the data sparse in the angular component. After 
applying such a transform the complete pressure data is recovered by solving a set of one 
dimensional .^^-minimization problems. This decomposition also makes our reconstruction 
algorithm numerically efficient. 

Although we focused on PAT using integrating line detectors, we emphasize that a similar 
framework can be developed for standard PAT based on the three dimensional wave equation 
and a two dimensional array of point-like detectors. In such a situation, finding a sparsifying 
transform is even simpler. Recalling that A^-shaped profile of the thermoacoustic pressure 
signal induced by the characteristic function of a ball suggests to use p as a sparsifying 

transform. 

Note that the recovery guarantees in this paper crucially depend on the choice of an appropri¬ 
ate sparsifying transform. In Subsection Kl.dl we proposed a candidate for such a transforma¬ 
tion that works well in our numerical examples. A more theoretical study of such sparsifying 
transforms and the resulting recovery guarantees for simple (piecewise constant) phantoms 
is postponed to further research. In this context, we will also investigate the use of different 
sparsifying temporal transforms, such as the ID wavelet transform in the temporal direction. 
Further research includes using a fixed number of detectors in each measurement process. 
This requires novel results for right fc-regular expander graphs and compressive sampling. 
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